#1.ANOVA检验
##数据处理##
library(openxlsx)
ADdata <- read.xlsx("E:\\class3\\13-刘睿佳-结果.xlsx")
rownames(ADdata)=ADdata[,1]   #将第一列作为行名
ADdata=ADdata[,-1]            #将第一列删除
ADdata[ADdata == 0]=NA        #将值为0的数据置为空
ADdatan <- apply(ADdata, 1:2, function(x) log2(x))   #对数化
colnames(ADdatan) <- gsub("\\d","",colnames(ADdatan))
group <- colnames(ADdatan)   #把样本信息分为ctl,asym,ad三组
table(group)

##ANOVA方差分析##
pvalue <- vector(mode="numeric",length=0)    #定义一个数值类型的空向量pvalue
for (i in (1:100)) {
  pro <- ADdatan[i,]        #提取第i行的蛋白质信息
  posctl <- group=='ctl'
  ctl <- pro[posctl]         #提取ctl组的数据
  ctlsum=sum(!is.na(ctl))    #计算ctl组非0个数
  posasym <- group=='asym'
  asym <- pro[posasym]         #提取asym组的数据
  asymsum=sum(!is.na(asym))    #计算asym组非0个数
  posad <- group=='ad'
  ad <- pro[posad]             #提取ad组的数据
  adsum=sum(!is.na(ad))        #计算ad组非0个数
  if(ctlsum>3 & asymsum>3 & adsum>3)  #对各组数据至少3个样本有数据的蛋白质进行ANOVA检验
  {result=oneway.test(pro~group)       
  pvalue[i]=result$p.value} else pvalue[i]=NA   #否则pvalue为空
}
pvalue
pofpro <- data.frame(pvalue,row.names = rownames(ADdata))
save(pofpro,file = "E:\\class6\\pvalue.RData")



#2.GO富集分析
library("clusterProfiler")
library(org.Hs.eg.db)
load("E:\\class5\\class5_volcano.RData")
data <- prostat[prostat$P<0.05,]     
allID <- data$ID[!is.na(data$P)]     #提取P<0.05的蛋白的名称
##将蛋白的名称数据进行格式转化##
ids<- bitr(allID, 
           fromType = "SYMBOL", #输入的格式为SYMBOL
           toType = c("GENENAME","ENTREZID"), #转化为ENTREZID格式
           OrgDb = 'org.Hs.eg.db')

##富集分析并绘制气泡图##
go.mf<- enrichGO(ids$ENTREZID,    #使用转换后的ENTREZID
                 OrgDb = org.Hs.eg.db,keyType ='ENTREZID', ont = "MF",   
                 pAdjustMethod = 'BH',pvalueCutoff = 0.05,
                 readable = T)           #分子功能（MF）富集分析
pdf("E:\\class6\\MF气泡图.pdf")    #存为pdf文件
dotplot(go.mf,showCategory = 20)   #绘制气泡图
dev.off()

go.cc<- enrichGO(ids$ENTREZID,    #使用转换后的ENTREZID
                 OrgDb = org.Hs.eg.db,keyType ='ENTREZID', ont = "CC",   
                 pAdjustMethod = 'BH',pvalueCutoff = 0.05,
                 readable = T)           #细胞组分（CC）富集分析
pdf("E:\\class6\\CC气泡图.pdf")    #存为pdf文件
dotplot(go.cc,showCategory = 20)   #绘制气泡图
dev.off()

go.bp<- enrichGO(ids$ENTREZID,    #使用转换后的ENTREZID
                 OrgDb = org.Hs.eg.db,keyType ='ENTREZID', ont = "BP",   
                 pAdjustMethod = 'BH',pvalueCutoff = 0.05,
                 readable = T)           #生物过程（BP）富集分析
pdf("E:\\class6\\BP气泡图.pdf")    #存为pdf文件
dotplot(go.bp,showCategory = 20)   #绘制气泡图
dev.off()